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We explore theoretically the physics of dynamic hysteresis for driven-dissipative nonlinear pho¬ 
tonic resonators. In the regime where the semiclassical mean-field theory predicts bistability, the 
exact steady-state density matrix is known to be unique, being a statistical mixture of two states: 
in particular, no static hysteresis cycle of the excited population occurs as a function of the driving 
intensity. Here, we predict that in the quantum regime a dynamic hysteresis with a rich phenomenol¬ 
ogy does appear when sweeping the driving amplitude in a hnite time. The hysteresis area as a 
function of the sweep time reveals a double power-law decay, with a behavior qualitatively different 
from the mean-field predictions. The dynamic hysteresis power-law in the slow sweep limit defines a 
characteristic time, which depends dramatically on the size of the nonlinearity and on the frequency 
detuning between the driving and the resonator. In the strong nonlinearity regime, the characteristic 
time oscillates as a function of the intrinsic system parameters due to multiphotonic resonances. We 
show that the dynamic hysteresis for the considered class of driven-dissipative systems is due to a 
non-adiabatic response region with connections to the Kibble-Zurek mechanism for quenched phase 
transitions. We also consider the case of two coupled driven-dissipative nonlinear resonators, show¬ 
ing that dynamic hysteresis and power-law behavior occur also in presence of correlations between 
resonators. Our theoretical predictions can be explored in a broad variety of physical systems, e.g., 
circuit QED superconducting resonators and semiconductor optical microcavities. 


I. INTRODUCTION 

Since the first experimental realization [l|, optical 
bistability has been the subject of many investigations. 
One of the major theoretical breakthroughs was accom¬ 
plished by Drummond and Walls, who provided an ex¬ 
act quantum solution for the steady-state density matrix 
of a sin^ mode driven-dissipative nonlinear optical res¬ 
onator y. One of the main conclusions is somewhat 
surprising at first sight since it reveals a unique steady- 
state solution while the semiclassical mean-field theory 
3 predicts a bistable regime with two stable branches. 
A large number of experimental studies in various sys¬ 
tems have shown optical hysteresis wcles (to give just 
a few recent references, see, e.g., 0 t 3) which seems in 
accordance with the semiclassical mean-field approach. 

This apparent contradiction can be reconciled by re¬ 
alizing that fluctuations (quantum or classical) induce 
switching between the two branches resulting in a unique 
stationary mixed state solution [Tol - [l^ . The switching 
is typically quantified by the lifetimes of the system in 
the different branches Without fluctuations, these 
lifetimes are infinite and produce the bistability at mean- 
field level. Note that the switching between the branches 
can be nicely visualized by considering individual quan¬ 
tum trajectories of the system [isl . [Tg. A similar be¬ 
havior is obtained if classical fluctuations are considered 
such as for example thermal fluctuations [T^ or a noisy 
drive d, [l^ • In general, it is possible to define a tran¬ 
sition point where the lifetimes of the two branches are 
equal. When the system is not at such transition point, 
one of the branches becomes increasingly more unsta¬ 
ble than the other. These lifetimes can however become 


extremely large with respect to all other timescales re¬ 
sulting in quasi-bistability (the solutions are metastable) 
and explaining the success of the mean-field theory to 
describe the hysteresis. In recent years, photonic res¬ 
onators with enhanced quantum nonlinearities have been 
developed [l^ , in particular using superconducting quan¬ 
tum circuits or semiconductor nanostructures as nonlin¬ 
ear media, paving the way to the experimental study of 
optical bistability in the quantum regime. 

In this paper, we report the surprising behavior of the 
hysteresis cycle in the quantum regime where the role 
played by quantum fluctuations and correlations becomes 
crucial. To explore such a physical problem, we have 
solved the time-dependent master equation for driven- 
dissipative nonlinear quantum resonators. By sweeping 
the drive amplitude in a finite time, we show that a dy¬ 
namic hysteresis is obtained, even when the steady-state 
quantum solution is unique. Within a time-dependent 
mean-field approach, the area of the hysteresis cycle con¬ 
verges to the finite steady-state mean-field value for in¬ 
finitely slow sweep with a deviation tending to zero as a 
power law, as it has been shown analytically and numer¬ 
ically for various classical models (20l - l^ . Here we show 
that the behavior of the hysteresis area in the quantum 
regime is qualitatively different in two ways: i) due to 
the uniqueness of the quantum steady-state solution, the 
area goes to zero in the limit of a very slow sweep; ii) the 
hysteresis area decays with increasing sweep time follow¬ 
ing a power law with an exponent that is different from 
the mean-field case for the same system parameters. We 
determine a characteristic time associated to the power- 
law decay and show its dramatic dependence on the size 
of the nonlinearity and frequency detuning of the driving. 
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In the regime of strong photon nonlinearity, an oscillating 
behavior of the characteristic time as a function of the 
system parameters is shown to be due to the quantization 
of the photon field. We show that the power-law behavior 
of the dynamic hysteresis can be captured analytically by 
determining a non-adiabatic response region. Concerning 
the experimental implementations, the presented physics 
is expected to be applicable to a broad range of models 
of driven-dissipative nonlinear photonic models that ex¬ 
hibit bistability at the mean field level. Examples are: 
the driven-dissipative Jaynes-Cummings model [^. 
an optomechanical cavity [ 2 ^ - 1 ^ . the driven-dissipative 
Dicke model [2^, and the micromaser [3l| . Recently 
there has been a strong increase of interest in coupled 
nonlinear photonic modes arranged in large lattice struc¬ 
tures for which the mean field approach typically pre¬ 
dicts bistability [H, . A better understanding of 

the role of quantum fluctuations is also crucial in the con¬ 
text of high speed optical switches of which the operation 
is based on optical bistability 


II. THEORETICAL FRAMEWORK AND 
PHYSICAL SYSTEMS 


Let us start by considering the quantum Hamilto¬ 
nian (/i = 1 ) for a single-mode boson field with boson- 
boson interaction and coherent driving (treated within 
the rotating-wave approximation): 


H{t) = Wcd^o -|- ^d^a’^aa -|- F{t)e -|- 

( 1 ) 

This is also the Hamiltonian of a single-mode cavity held 
with a dispersive (Kerr) optical nonlinearity. Here, the 
boson operator a (d^) annihilates (creates) an excitation 
in the cavity, while ujc is the photon mode frequency, 
U quantihes the photon-photon interaction, F(t) is the 
time-dependent driving amplitude, and ujp is the driving 
frequency. 

The dynamics with dissipation is described by the 
Lindblad master equation for the density matrix p(t): 


dKt) 

dt * 


p,H{t) -I- ^ (1 -I- Tith) (2dpd^ — a^ap — pdfd) 
+ 2 ''^th pa — aa}p — pddf) , (2) 


where the term proportional to the commutator 


p,H{t) 


describes the quantum dynamics due to the Hamiltonian. 
The quantity 7 is the dissipation rate due to the cou¬ 
pling to the environment, nth the mean-number of ther¬ 
mal excitations at the resonator frequency Wc, namely 
nth = — 1 ) ^, with /3 = T the bath 

temperature and ks the Boltzmann constant. From the 
density matrix the expectation value of an observable 


O can be calculated as (O) = Tr 


Op 


, where Tr de¬ 


notes the trace. For the calculations, we will work in the 


frame rotating at the pump frequency for which the de¬ 
tuning A = Wp — Wc is the relevant parameter. Notice 
that even in the rotating frame the Hamiltonian remains 
time-dependent due to time-dependence of the driving 
amplitude Fit). 

For comparison, the mean-field equation for the coher¬ 
ent field a{t) = (a) reads: 

=(^uz.-tl + U \af) a + F(t)e—(3) 

In the following, we will systematically compare the pre¬ 
dictions of the exact solutions of the master equation ([ 2 ]) 
to those obtained within the mean-field approximation 
in Eq. ([3|). 

The present theoretical model is rather general and can 
be obtained, e.g., in a system consisting of a coherently 
driven linear cavity coupled to an ensemble of two-level 
atoms in the dispersive limit (large detuning between 
cavity and atom resonance frequencies with respect to 
the coupling strength) [d^. Novel quantum optical sys¬ 
tems with large nonlinearities such as superconducting 
quantum circuits and semiconductor microcavities have 
emerged in recent years M- In semiconductor pillars 
including quantum wells, a normalized U /y up to a few 
percent is within present capabilities. For these systems, 
a temperature T = 4K (liquid helium bath) together with 
a cavity photon energy of 1.5 eV give [3ujc ^10^ and a 
negligible number of thermal excitations: nth ~ 0. In 
the context of circuit QED where a Josephson junction is 
used to introduce an effective nonlinearity [1, M 113) 113 > 
much larger nonlinearities can be be achieved with the 
possibility to reach the hardcore boson limit \U\/^ ^ 1 
[50l |. A typical dilution fridge temperature of 50 mK and 
a resonator frequency Uc/ (Stt) = 5 GHz corresponds to 
Puic — 4.8 and nth — 0.008. Hence, thermal effects are 
in general small, but not completely negligible in circuit 
QED. 

Note that here for simplicity we will not consider sys¬ 
tems with absorptive optical nonlinearity and bistabil¬ 
ity (i.e., a cavity photon mode resonant with an elec¬ 
tronic excitation like in the celebrated Jaynes-Cummings 
model) where quantum fluctuations also induce switching 
between the semiclassical branches [H, [13 . 

In order to solve numerically the master equation ([2) , 
we have expressed the time-dependent density matrix in 
the basis of Fock number states \n), namely Pn,m{t) = 
(n |/5(t)| m). Convergence of the results have been care¬ 
fully checked by increasing the cut-off number of photons. 
With this numerically exact integration method, we can 
typically explore regimes with number of photons up to 
few tens. 


HI. NUMERICAL RESULTS 

In order to study dynamical hysteresis phenomena, we 
consider a triangular modulation of the drive amplitude. 
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FIG. 1. (a) The photon population n and (b) the 

second-order correlation function versus the driving ampli¬ 
tude F (units of 7 ) for a single-mode driven-dissipative quan¬ 
tum resonator with a nonlinearity U = O.I 7 and detuning 
A = 27 . In panel (a), the steady-state mean-field (MF) re¬ 
sult and the quantum steady-state solution (SS) from Ref. 
Q| are presented. The other two curves are dynamic hys¬ 
teresis cycles predicted by the time-dependent quantum mas¬ 
ter equation obtained by using two different sweep times ts 
{ts/AF = 10 / 7 ^ for the curve with the largest hysteresis cy¬ 
cle and ts/AF = 20 / 7 ^ for the smaller one). In panel (b) 
the steady-state solution is shown together with the result 
for a time-dependent sweep with ts/AF = lO/y^ (the arrows 
indicate the direction of the sweep). 


namely consisting of one sweep from Fq to Fq -f AF and 
one from Fq + AF back to Fq : 


f t — 9t 

F{t) = Fo + -AF9{ts -t)- -^AFe{t - F), (4) 


where 9{t) is the Heaviside step function and the time pa¬ 
rameter ts is the sweep time. In practice the parameters 
Fq and AF are always chosen such that the sweep covers 
the full range of the hysteresis. The master equation is 
solved in the time interval from t = 0 to t = 2ts with 
the steady-state solution at the pump intensity Fq as an 
initial condition. Note that the presented results are in 
the zero temperature limit (/3wc —^ +00 and nth -9- 0 ) 
unless explicitly stated otherwise. 



FIG. 2. (a) The area A of the hysteresis loop as a func¬ 

tion of the sweep time ts (units of AF/ 7 ^) for different tem¬ 
peratures (from bottom to top the thermal population nth 
is 0.2,0.1,0.05 and 0, corresponding respectively to Pujc — 
1.8, 2.4, 3 and + 00 ), together with the result from the mean- 
field approximation (MF) for U = 7/2 and A = 27 . The 
full lines are power law fits to the different limiting regimes 
for which two separate power laws are observed. For large ts 
we find the behavior A (x t~^ while for small values of ts we 
find: A oc tj** with a coefficient b that depends on the sys¬ 
tem parameters. For the mean-field result we find an overall 
good agreement with {A — Aq) oc tj^^^ with Ho > 0 the static 
hysteresis area. The characteristic timescale r, as determined 
from the behavior A — (ts/(rAF))“^ for large ts, is shown in 
(b) as a function of the nonlinearity U (units of 7 ) for differ¬ 
ent values of the detuning A and in (c) as a function of the 
detuning A (units of 7 ) for different values of the nonlinear¬ 
ity U. Note the oscillating behavior with minima satisfying 
the n-photon resonance conditions: Un{n — l)/2 = nA. The 
characteristic timescale r in Fig. (a) is II 5/7 . 


A. One resonator 


In Fig. [T] (a), we consider results for the excited popu¬ 
lation n = (a/d) using parameters for which the steady- 
state semiclassical mean-field (MF) solution exhibits the 
well-known bistability of the coherent population jap as 
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a function of the driving amplitude (three branches with 
the middle one unstable). In stark contrast, the steady- 
state (SS) solution [2] of the master equation shows no 
hysteresis cycle for n. However, if we consider a time- 
dependent solution of the master equation with finite 
sweep time tg, a clear dynamic hysteresis is found for 
n{F{t)). The area of the dynamic hysteresis loop ob¬ 
tained with the solution of the master equation decreases 
for increasing tg (slower sweep). In the adiabatic limit of 
an infinitely slow sweep {tg —>■ -|-oo) the hysteresis disap¬ 
pears and the steady-state density-matrix solution of the 
master equation is recovered. In Fig. [1] (b) the normal¬ 
ized second-order correlation function = (d'^d^dd) / 

is presented as a function of the driving amplitude to¬ 
gether with the steady-state result from Ref. [^. At 
a transition point, a sharp peak with significantly 
larger than 1 occurs, a feature that cannot be captured 
by mean-held theory for which g^^^ = 1. For the time- 
dependent solution of the master equation, we hnd that 
this peak is shifted with respect to the steady-state so¬ 
lution to the region where the transition is seen in the 
density and it is more (less) pronounced for decreasing 
(increasing) F. 

In order to study quantitatively the properties of dy¬ 
namical hysteresis, we will evaluate the area A of the hys¬ 
teresis loop for a sweep with population n-\-{F) obtained 
for increasing driving amplitude F and population ni{F) 
achieved for decreasing F: 

/*Fo+AF 

A= dF\n^{F)-n^iF)\. (5) 

JFo 

In Fig. [2] (a) the area A is plotted as a function of the 
sweep time tg for different temperatures together with 
the result predicted by the mean-field equation ([3]). Our 
results show that for relatively fast sweeps (small tg) a 
reasonable agreement is found between the exact solu¬ 
tion and the time-dependent mean-field result while for 
slower sweeps qualitatively different scaling laws are ob¬ 
served (the full lines in Fig. [2] (a) are power-law fit¬ 
ting curves). In the regime where the mean-field steady- 
state solutions exhibit bistability (A > the time- 

dependent mean-field equation predicts a dynamic 
hysteresis area which can be well fitted by the expression 
(A — Ao )/7 oc with Aq the hysteresis area in the 

adiabatic limit {tg —^ -l-oo). This is in agreement with 
the analytic result from Ref. for a similar mean- 

field equation. In stark contrast, we find that the exact 
solution of the master equation has a double power-law 
behavior: for large tg the area scales as A oc tj^. In Sec¬ 
tion [IVl we report an analytical derivation of the power 
law based on the determination of a non-adiabatic region. 
In Appendix |A] , we present additional numerical results 
using a quasi-adiabatic approximation. In the regime 
where mean-field predicts no bistability (A < \/i/2'y) 
the dynamic mean-field result also exhibits the power law 
A oc for slow sweeps. Therefore, in contrast to the 
exact quantum solutions, we emphasize that the mean- 
field approach predicts different scaling laws depending 


on the frequency detuning, a result in agreement with 
Refs. [U, treating mean-field models. Furthermore, 
we point out that our time-dependent quantum result 
exhibits another power law at small values of tg with a 
coefficient that depends on the system parameters (see 
Fig. [21(a)). Note that this kind of double scaling law for 
the hysteresis area has also been observed in the context 
of dynamic transitions with magnetic materials [5l|, . 

Moreover, we see that the presence of a moderate thermal 
population (typical values of circuit QED experiments) 
gives the same power law behavior and just a moderate 
decrease of the hysteresis area as the temperature is in¬ 
creased. This is expected since thermal fluctuations also 
contribute to switching between the two branches. 

The power-law behavior allows us to determine a char¬ 
acteristic timescale r on which the hysteresis size will 
change significantly, namely through the expression A = 
(ts/(TAF))“^, in the regime of large tg (see Fig. |2] (a)). 
For a sweep with tgj/AF ^ r the quantum fluctuations 
are expected to induce a significant deviation from the 
mean field result (as can also be seen in Fig. [3] (a) where 
t/j = 115). The characteristic timescale r is presented 
as a function of the nonlinearity for different values of 
the detuning in Fig. |2] (b) and as a function of the de¬ 
tuning for different values of the nonlinearity in Fig. |2| 
(c). We point out that the characteristic time r can be 
orders of magnitude larger than the resonator lifetime 
1/7 for small nonlinearities and/or large detuning. As 
a function of the detuning, an overall exponential in¬ 
crease of the characteristic hysteresis time is observed. 
Note that the dynamic hysteresis survives in the hard¬ 
core limit ([/ —>- 00 ) with r converging to a finite value. 
Furthermore, a superimposed oscillating behavior of the 
characteristic time is predicted, which becomes more pro¬ 
nounced as the detuning or nonlinearity is increased. The 
minima correspond to system parameters satisfying the 
resonance condition Un{n — l)/2 = nA, with n a posi¬ 
tive integer (giving the sequence U = 2A, A, 2/3A,...). 
These are the n-photon resonances [s^, obtained when 
the energy of n pump photons is equal to the energy of 
n interacting photons in the resonator. 

B. Two coupled resonators 

Now, we consider the case of two identical resonators 
coupled by the following hopping Hamiltonian: 

Fjiop — J ^a^a2 F h.c^ , (6) 

with J the hopping parameter. Note that this model 
is currently experimentally realized, e.g., with semicon¬ 
ductor microcavities [H, and with circuit QED res¬ 
onators [ 5 ^. For this system, for sake of simplicity, we 
have considered the case when each resonator has the 
same driving term. Since adding a second resonator 
squares the dimension of the Hilbert space, in order to 
have exact results with arbitrary precision, we are re- 
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FIG. 3. Results for a system consisting of two coupled iden¬ 
tical resonators, (a) The inter-resonator correlation function 
as a function of the pump amplitude F (in units of 7 ) 
during a sweep of the driving amplitude. Results are shown 
for the steady state (SS) and for a dynamic hysteresis with 
ts/(AF’) 7 ^ = 30 (the arrows indicate the direction of the 
sweep), (b) The hysteresis area A associated to the photon 
population in one resonator as a function of the sweep time to¬ 
gether with the mean-field (MF) result. The full lines are fits 
obtained with the following scaling laws: A oc for large ts 
and {A — Ao)/j oc for the mean-field result. The other 
system parameters A, U and J are specified in each panel. 


stricted to lower photon numbers with respect to the sin¬ 
gle resonator case. A negative detuning A is considered 
which is blue-detuned with respect to the linear reso¬ 
nance corresponding to the ’bonding’ single-particle state 
(having eigenfrequency lu = uic — J)- As in the single¬ 
resonator case, we find a dynamic transition around the 
regime where mean-field predicts bistability with the hys¬ 
teresis cycle and a peak in the local second-order corre¬ 
lation function at the transition (qualitatively the same 
as observed for a single cavity in Fig. [T|). We also exam¬ 
ined the inter-cavity normalized second-order correlation 
function: g ^2 = /{nin 2 )^ which is presented 

in Fig. [3] (a) for a temporal sweep and for the steady- 
state. Also for this non-local correlation function a peak 
is observed at the transition. In Fig. |3] (b) we have pre¬ 
sented the hysteresis area, as defined in Eq. ( 0 , where 
the population of one resonator mode is used (due to the 
symmetry the populations are equal in both resonators). 
This reveals qualitatively the same power laws as for a 
single cavity: A oc for large ts and [A — Aq) oc 
for the dynamic mean-field result. These results for two 
coupled cavities are relevant considering the recent in¬ 
terest in strongly correlated photonic phases in arrays 


of cavities (see for example [l^ Isl - fsM I40l - l4l |). It was 
shown that applying the Gutzwiller mean-field approach 
to a lattice of nonlinear driven-dissipative cavities pre¬ 
dicts bistability [13,113 ■ Our exact results for the two- 
resonator system show that while the exact solution has 
no static hysteresis, a dynamic hysteresis does emerge 
in the quantum regime. A natural next step to be pur¬ 
sued in the future is a study of the dynamic hysteresis 
for larger arrays of coupled resonators. In particular, the 
dependence of the dynamic hysteresis on the inevitable 
presence of disorder, the effect of multimode dynamics 
and quasi-continuous spectrum of states is an open prob¬ 
lem that would be interesting to explore. 

IV. ANALYTICAL SCALING BEHAVIOR AND 
CONNECTION WITH KIBBLE-ZUREK 
MECHANISM 

In the previous section, we have presented a compre¬ 
hensive set of numerical solutions of the master equation 
showing the rich properties of dynamic hysteresis for a 
driven-dissipative nonlinear quantum resonator. In this 
section, we present an analytical demonstration of the 
power-law behavior and of the exponent. Qualitatively, 
we show that the dynamical hysteresis is due to a non- 
adiabatic response of the considered system when the 
driving field is swept around the bistability region. Note 
that the approach presented here can be applied to a 
generic driven-dissipative system. 

When changing in time one parameter of an Hamil¬ 
tonian system, by definition the response becomes non- 
adiabatic when the time scale of the change is much 
shorter than the time scale of the system internal dy¬ 
namics. Such a time is proportional to the inverse of 
the energy gap between the ground state and the excited 
state manyfold. In the case of quantum phase transitions, 
the energy gap vanishes at the critical point (softening of 
the excitation mode) leading to a divergence of the cor¬ 
responding internal dynamics time scale (critical slowing 
down). Therefore, when crossing a critical point, there is 
always a non-adiabatic response region around the tran¬ 
sition. This property is at the heart of the celebrated 
Kibble-Zurek mechanism for the formation of topologi¬ 
cal defects in quenched quantum phase transitions [56l - 
Issjl (see for example Ref. for a review). Recently, 
the Kibble-Zurek mechanism has been examined for the 
non-dissipative quantum Rabi model which consists of a 
zero-dimensional photonic mode coupled to a single two- 
level system [^ . 

Since the class of systems we are studying in the 
present paper is of the driven-dissipative kind, the gap of 
the Hamiltonian is not at all the relevant quantity. What 
is relevant here is the spectrum of the Liouvillian super¬ 
operator L associated to the master equation dtp = Lp. 
We consider the eigenvalue equation for the Liouvillian 
superoperator: 

Lp\ = Xpx, 


(7) 
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where the eigenvalue A is in general complex. The steady- 
state density matrix corresponds to the eigenvalue A = 0. 
Note that the real part of A represents the damping of 
the excitation mode, while the imaginary part represents 
the oscillation frequency with respect to the steady-state. 
Here, we will focus on the eigenvalue with the smallest 
non-zero real part since it is the least damped and de¬ 
termines the asymptotic relaxation to the steady-state. 
In Fig. |4] the real and imaginary part of this eigenvalue 
are presented as a function of the drive amplitude for 
U = O.I 7 and A = 27 . Note that the the frequency 
gap (imaginary part) is 0 in a finite interval around the 
value Fc where the damping (real part) is also strongly 
suppressed, albeit reaching a finite minimum. This re¬ 
veals that such a mode is soft and diffusive, i.e. because 
it is degenerate in frequency with the steady-state and 
has a finite damping. This kind of soft diffusive modes 
can appear in driven-dissipative systems (in the case of 
Bogoliubov excitations of driven-dissipative systems, see 
[ 1 ^ 1 . Hence, for a sweep of F around Fc the response 
of the system is expected to have a non-adiabatic con¬ 
tribution because of such frequency degeneracy. From 
this diffusive soft mode we determine the relaxation time 
tr = —l/Re[X\. Note that the so-called tunneling time 
tt of bistability [lol - [l^ corresponds to the maximal value 
of T/j, at the transition point. 




FIG. 4. The real (a) and the imaginary (b) part of the Liou- 
vilian eigenvalue A (in units of 7 ), corresponding, respectively, 
to the damping rate and the frequency of the excitation mode. 
In particular, we consider the least damped mode (different 
from the steady-state corresponding to A = 0) as a function 
of the drive amplitude F (in units of 7 ) for U = O.I 7 and 
A = 27 . Around the transition point (at Fc ~ 87 ) the damp¬ 
ing rate (real part) is strongly suppressed, while the imaginary 
part is exactly zero, indicating the presence of a soft diffusive 
mode. Away from the transition region there are two sym¬ 
metric least damped modes with equal damping rates but 
opposite frequencies (the imaginary parts). 


For a time-dependent sweep of the drive amplitude, we 
introduce the distance e(t) from the value Fc, namely: 

e{t) = Fc-F{t). (8) 

We consider now a sweep of F{t) linear in time from 
Fc — AF/2 to Fc + AF/2 with total time duration tg. 


FIG. 5. (a) The relaxation time tr (the curve peaked around 
Fc « 87 ) is presented as a function of the drive amplitude 
together with the sweep timescale Tg for two sweeps with dif¬ 
ferent speeds for U — O.I 7 and A = 27 . The non-adiabatic 
region with width 5F around Fc is indicated for the fastest 
sweep, (b) The width of the sweep SF as a function of the 
sweep time together with the two power laws, (c) The tunnel¬ 
ing time Tt (full line) and the characteristic time r (dashed 
line) versus the detuning A for nonlinearities [ 7/7 = 4 (lower 
curves) and [ 7/7 = 1 (upper curves). 


The normalized sweep rate reads (s^ 

e{t) tg \Fc-F{t)\ T,’ 

which defines a sweep time scale Tg{F) = tg\Fc — F\/AF, 
which is plotted in solid thin line in Fig. [S] (a). Note 
that an alternative derivation of this expression can be 
obtained by equating the time from the transition point 
Fc to the instantaneous relaxation time Note that 
the sweep time scale Tg differs from the duration tg, be¬ 
cause it contains also information on how much the driv¬ 
ing amplitude is changed with respect to the transition 
point. 

By generalizing the criterion used for the Kibble-Zurek 
mechanism in the context of equilibrium phase transi¬ 
tions, we expect that our driven-dissipative system enters 
a non-adiabatic regime when the sweep timescale Tg be¬ 
comes smaller than the relaxation time tr. In Fig. [5] (a), 
we plot the relaxation time tr and the sweep timescale Tg 
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as a function of the driving amplitude F for two different 
sweeps {tgj^/AF = 10^ and 10^). The corresponding 
non-adiabatic region of width 6F around Fc is indicated 
in Fig. [5] (a) for the fastest sweep {tsj^/AF = 10^). In 
Fig. [5] (b) the width SF of the non-adiabatic region is 
presented versus the sweep time duration ts, showing a 
double power-law, as also found for the area of the dy¬ 
namic hysteresis. Furthermore, for slow sweeps the decay 
is proportional to the same exponent as found be¬ 
fore. In the non-adiabatic region the system does not 
have the time to relax to the steady-state, resulting in 
an hysteretic behavior. The area of the hysteresis loop 
is therefore linked to the width SF of the non-adiabatic 
region, as confirmed by our numerical results. 

For ts —>■ +00 (slow sweep limit), the non-adiabatic 
region size SF —> 0. Hence, in this slow sweep limit, 
Tfi ^ tt, the maximum of the relaxation time. From Eq. 
(|ni) and imposing Tg = tb, ~ tt, we get the asymptotic 
expression for dF, namely: 

SF = 2TT(ts/AF)-\ (10) 

This formula is in excellent agreement with the results 
in Fig. [S] (b). We also note that for smaller values of F, 
this approximation fails and so a different power-law is 
expected, depending on the shape of the relaxation time 
vs F. This shows that the exponent —1 of the slow sweep 
power law occurs because the real part of the Liouvil- 
lian eigenvalue remains finite, but much smaller than all 
other characteristic frequencies of the system. The power 
law exponent is thus expected to change when also the 
real part of the Liouvillian eigenvalue vanishes; i.e. when 
tt —» oo, corresponding to a dissipative phase transition 
|6ll |. This might occur in the thermodynamic limit of a 
large lattice of coupled driven-dissipative resonators. In 
Fig. [5] (c) the tunnelling time tt is compared with the 
characteristic time r (see previous section and Fig. [Hi as 
a function of the detuning A for two values of the nonlin¬ 
earity. This reveals qualitatively similar behavior with an 
overall exponential increase as a function of the detuning 
and oscillations due to the multi-photonic resonances. 

For conservative systems with a finite energy gap the 
Kibble-Zurek mechanism breaks down for slow sweeps 
since the evolution becomes adiabatic [^, In 

this case an effective description is provided by the 
Landau-Zener approximation for the evolution of a sys¬ 
tem through an avoided energy crossing [gj, • Apply¬ 
ing the Kibble-Zurek mechanism results in a good agree¬ 
ment with the Landau-Zener result only for sufficiently 
fast sweeps [g^, [g3- Note that the Landau-Zener for¬ 
mula for a dissipative excited state does not depend on 
the decay rate (gg| and its applicability is connected to 
the existence of a finite gap for the frequency. 

For the considered dissipative system on the other 
hand we find that the scaling laws based on a Kibble- 
Zurek-like approach for the non-adiabatic regime agree 
with the numerical results for the hysteresis area, also 
in the slow sweep limit. This shows that an adiabatic 
regime is never reached, no matter how slow is the sweep. 


At first sight this might seem in conflict with the results 
for the Kibble-Zurek mechanism for conservative systems 
since the real part of the Liouvillian gap remains finite. 
However, for dissipative systems it is the imaginary part 
of the Liouvillian eigenvalue that gives the excitation en¬ 
ergy. 


V. CONCLUSIONS 

We have investigated the time-dependent exact so¬ 
lutions of the quantum master equation for driven- 
dissipative nonlinear quantum resonators, thus includ¬ 
ing the role of quantum fluctuations and correlations. In 
particular, we have focussed on the regime where the 
semiclassical mean-field approximation predicts bistabil¬ 
ity and investigated temporal sweeps of the drive ampli¬ 
tude revealing dynamic hysteresis loops. The hysteresis 
behavior, typically attributed to the semiclassical mean- 
field approach, was found to survive in the regime of small 
photon numbers and strong quantum fluctuations. The 
time-dependent quantum solution, in contrast to predic¬ 
tions of mean-field approaches, shows that the hysteresis 
area as a function of the total sweep time has a double 
scaling law. These results have been shown to be robust 
with respect to thermal excitations for typical experi¬ 
mental temperatures. We have determined a character¬ 
istic time associated to the power-law decay of the dy¬ 
namic hysteresis area, showing a rich behavior as a func¬ 
tion of the nonlinearity and of the frequency detuning. 
We have also considered two coupled driven-dissipative 
resonators, demonstrating that dynamic hysteresis and 
power-law decay occur also in presence of inter-cavity 
correlations. Importantly, we have demonstrated that 
the dynamic hysteresis is associated to a non-adiabatic 
response region with connections to the Kibble-Zurek 
mechanism for quenched phase transitions. We have been 
able to describe analytically the power-law behavior with 
scaling arguments and shown the role of a soft diffusive 
mode, i.e. having zero excitation energy, but a finite 
damping. This is a general picture, which is expected to 
apply to a broad class of driven-dissipative quantum sys¬ 
tems. These exciting results can be a motivation to fur¬ 
ther investigate larger arrays of cavities at the dynamic 
transition. Given the emergence of several interesting 
systems with controllable quantum optical nonlinearities, 
the present predictions should stimulate exciting studies 
of dynamic hysteresis in the quantum regime. 
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Appendix A: Quasi-adiabatic approximation 

In this Appendix a quasi-adiabatic approximation is 
applied to study the dynamic hysteresis. From the Lind- 
blad master equation the following equation of mo¬ 
tion for the density n can be derived: 

f)'Yl 

— =-^n-i{F{t){a^) - F*{t){a)) . (Al) 

As in the main text, we consider a linear sweep of the 
driving amplitude F{t) (see Eq. (|1])). For an adiabatic 
sweep [ts —>■ oo), the time-derivative tends to zero and 
all expectation values converge to the steady-state (SS) 
solution: 


'^ts—^oc t U'5'5', (-^^) 

t {cl)sS' (-^^) 

An analytical expression for the steady-state correlation 
functions was derived in Ref. Q. The result for the 
boson coherence (a) is (for sake of clarity we write the 
dependence on the driving amplitude F explicitly): 


{a)ss{F) 


F F{l + c,c*,8 
A + *7/2 jr(i + c,c*,8 



(A4) 


with c = — 2(A -I- i^l2)lU and F{c,d,z) the hypergeo¬ 
metric function: 


-A(c, d,z) 


^ r(c)r(d) z" 
r(c -|- n)r(d -I- n) n! ’ 

n—0 


(A5) 


r being the gamma special function. If the sweep is per¬ 
formed sufficiently slowly, it is interesting to see what are 
the predictions of a quasi-adiabatic approximation. This 
corresponds to using the exact steady-state solution for 
(a): (a) ^ (d) Ss{F{t)), where the time-dependent driv¬ 
ing amplitude F{t) is used in Eq. (IA4I1 . By performing 
this substitution in the equation of motion for the den¬ 
sity mil and numerically solving the differential equa¬ 
tion, the dynamic hysteresis can be examined within the 
quasi-adiabatic approximation. In Fig. [5] the resulting 
hysteresis area is compared to the exact numerical result 
for a set of parameters. The results clearly deviate, re¬ 
vealing that the quasi-adiabatic approximation does not 
capture the full dynamics. However, the same power law 
behavior A oc is found for sufficiently slow sweeps. 
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